EDA workflow - #804

Description

This is a first step to make EDA workflow run within SWELL framework. These 32 ensemble members from @rtodling run independently on each node. Both GSI B and all observation type from 3dvar-atmos are used here. The single node run for each member is made possible with the current setting: C90 resolution using thinned observations (rejecting 80% of obs at random).
The memory usage on each node is roughly 72% (There are still details to work out and improve).

The overall workflow:

The EDA-3D implementation here is a jedi block-EDA version. It is a very mechanical version of the EDA method.
A conventional 3D-Var minimization is invoked for each member state and the resulting ensemble anlysis states are final products of 3D-EDA. The observation perturbation term is added to the cost function when the member is a perturbed state.

1. The driver for this mechnism is in the .cylc file:
   

   swell/src/swell/suites/eda/flow.cylc

   Lines 79 to 81 in 35e846c

   	{% for i in range( 1, models[model_component]['ensemble_num_members'] + 1 ) %}
   	SP2-{{model_component}} => RunJediEdaExecutable_mem{{i}}-{{model_component}} => SP3-{{model_component}}
   	{% endfor %}

   
   

   swell/src/swell/suites/eda/flow.cylc

   Lines 170 to 179 in 35e846c

   	{% for i in range( 1, models[model_component]['ensemble_num_members'] + 1 ) %}
   	[[RunJediEdaExecutable_mem{{i}}-{{model_component}}]]
   	script = "swell task RunJediEdaExecutable $config -d $datetime -m {{model_component}} -imem {{i}}"
   	platform = {{platform}}
   	execution time limit = {{scheduling["RunJediEdaExecutable"]["execution_time_limit"]}}
   	[[[directives]]]
   	{%- for key, value in scheduling["RunJediEdaExecutable"]["directives"][model_component].items() %}
   	--{{key}} = {{value}}
   	{%- endfor %}
   	{% endfor %}

Each job submission completes one cost function minization for one member:
"swell task RunJediEdaExecutable $config -d $datetime -m {{model_component}} -imem {{i}}"

2. The observation thinning can be disabled or tuned by setting obs_thinning_rej_fraction from zero (no thinning) to 1 (reject all).
   

   swell/src/swell/suites/eda/suite_config.py

   Line 107 in 35e846c

   qd.obs_thinning_rej_fraction(0.8),

Thinning can also be turned off by commenting out this line:

3. Each EDA member run uses its own copy of observation files to prevent HDF5 file lock when multiple jobs are accessing a single file. Hence cleaning up the directory for each 3D-Var run is important:
   

   swell/src/swell/suites/eda/flow.cylc

   Lines 88 to 89 in 35e846c

   	# CleanupEda directory
   	SP4-{{model_component}} => CleanEdaFiles-{{model_component}}

The diagnostic tools:

For diagnostic purpose, we have designed handy codelets to compute mean/variance and diffstate (for increment, e.g.) in swell.

A. To compute the mean and variance of a set of states including bkg and analysis,
user can define the file specs for a list containing multiple entries where each can contain
its own input file, output file, and grid_type.

There are two ways to create the configure file.

Method-1. modify the EDA suite config file:

Method-2. create a section in override.yaml

models:
  geos_atmosphere:
    ensmeanvariance_spec:
    - state: bkg
      fn_input: ebkg/mem%mem%/geos.mem%mem%.%yyyy%mm%dd_%hh%MM%ssz.nc4
      fn_output_mean: geos.prior.mean
      fn_output_variance: geos.prior.variance
      grid_type: [ 'latlon', 'cs' ]
    - state: analysis
      fn_input: analysis/mem%mem%/eda.ana.mem%mem%.%yyyy%mm%dd_%hh%MM%ssz.nc4
      fn_output_mean: eda.ana.mean
      fn_output_variance: eda.ana.variance
      grid_type: [ 'latlon', 'cs' ]

and run: swell create eda_atmos -o override.yaml to override the config default

Finally to compute the mean and variance, run
swell task RunJediEnsembleMeanVariance PATHTO/experiment.yaml -d $date -m geos_atmosphere

B. To compute the difference between two states, similar approaches can be used.

Method-1. modify the EDA suite config file:

Method-2. create an override section in override.yaml

models:
  geos_atmosphere:
    ensmeanvariance_spec:
    diffstates_spec:
      state1:
        fn_input: geos.prior.mean.%yyyy%mm%dd_%hh%MM%ssz.nc4
      state2:
        fn_input: eda.ana.mean.%yyyy%mm%dd_%hh%MM%ssz.nc4
      state_diff:
        fn_output: eda.mean-inc
        grid_type: [ 'cs', 'latlon' ]  
      state_type: ensemble

and run: swell create eda_atmos -o override.yaml to override the config default

Finally to compute the diff states, run
swell task RunJediDiffstates PATHTO/experiment.yaml -d $date -m geos_atmosphere

How to run this block EDA code:

swell create eda_atmos -o override.yaml  -p nccs_discover_sles15 -k
swell launch   your_eda_suite
[override.yaml contains experiment_root and experiment_id]

Example output: T mean increment at 500 hPa :

The work is advised by @rtodling.
@mranst makes the flow.cylc loop work for ensemble members.

Timing results

In progress

Additional features:

1. In analogy to fv3jedi ctest for 3d-EDA, we have implemented a choice to use hybrid-B matrix in 3D-Var. Since this method is not tested, we have coded this function as an option for future use.

https://github.com/GEOS-ESM/swell/blob/35e846ce2371c9ac360ed13cf852082928bf922e/src/swell/configuration/jedi/interfaces/geos_atmosphere/model/background_error_eda_hybridB.py

A technical interesting part of the code is that I am using a string template in .py file, which is very useful here because it avoids writing 72 lines of pure real numbers into python code, and because we will need to tune and update this set of coefficient quite often later on, presumably. This is viewed as a transition back from python code to yaml style.

Issues

1. This code skips obs perturbation when the obs space contains covariance for R, because we did not find the corresponding source code in JEDI for this function.

2. People should note that:
   Adding air_pressure_levels to the list of analysis variable will break the 3D-Var code with hybrid-B. (I assume there is an issue with bump code) So I did not include it in the interface model config setup for EDA. Here my assumption is that always having air_pressure_at_surface included in analysis variable can be sufficient.

To do list:

1. We can eliminate hofx output, increment output in the 3D-Var run for each member, but we want to compare with output from more advanced EDA method.